“Baby Mental Life: Study 2” was conducted on MTurk on 2018-08-04.
Our planned sample was 300 participants, and we anticipated that roughly 80% of recruited participants would pass all of our attention checks, so we initially recruited 378 participants (on the idea that ~80% of 378 ~ 300 participants; note that for administrative purposes we need to recuit participants in batches that were divisible by 9). After filtering out participants who failed at least one of our attention checks, we ended up retaining fewer than 300 participants, so we recruited an additional 16 participants for a total of 394 people recruited. At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.
In the end, we ended up with a sample of 304 participants who passed our attention checks, 237 of whom came from unique GPS coordinates.
For this first pass, these data exclude participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.
Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 20 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable).
For more details about the study, see our preregistration here.
Here we run some exploratory analyses on these data.
# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
# set theme for ggplots
theme_set(theme_bw())
# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/table_fun.R")
source("./scripts/data_prep.R")
NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"
First, I’ll do the factor analysis, and get a list of the top 5 items by factor:
# load in S1 efa in case we need it
efa_S1 <- readRDS("../study 1/s1_efa.rds")
demo_S1 <- read.csv("../study 1/s1_demo.csv")
# conduct S2 efa
efa_S2 <- fa(d_all, nfactors = 4, rotate = "oblimin", fm = "minres",
scores = "tenBerge", impute = "median")
Loading required namespace: GPArotation
table_fun(efa_S2, num_items = 5, pos_abs = "abs")
| capacity | loading |
|---|---|
| Factor 1 | |
| feeling_distressed | 0.82 |
| feeling_overwhelmed | 0.82 |
| feeling_frustrated | 0.78 |
| feeling_helpless | 0.76 |
| feeling_lonely | 0.60 |
| Factor 2 | |
| having_self_control | 0.96 |
| controlling_their_emotions | 0.93 |
| telling_right_from_wrong | 0.93 |
| planning | 0.89 |
| reasoning_about_things | 0.89 |
| Factor 3 | |
| getting_hungry | 0.83 |
| feeling_pain | 0.79 |
| feeling_tired | 0.67 |
| hearing_sounds | 0.58 |
| feeling_physically_uncomfortable | 0.49 |
| Factor 4 | |
| feeling_happy | 0.82 |
| finding_something_funny | 0.81 |
| feeling_excited | 0.79 |
| loving_somebody | 0.64 |
| learning_from_other_people | 0.51 |
Now I’ll use this to define categories of mental capacities:
# get items by factor
factors_S2 <- efa_S2$loadings[] %>%
data.frame() %>%
rownames_to_column("capacity") %>%
gather(factor, loading, -capacity) %>%
group_by(capacity) %>%
top_n(1, loading) %>%
ungroup() %>%
select(-loading) %>%
mutate(factor_names = recode_factor(factor,
"MR1" = "Negative emotions",
"MR2" = "Cognition & control",
"MR3" = "Bodily sensations",
"MR4" = "Positive/social emotions"),
factor = factor(factor))
# make new dataframe
d_cat <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(d_demo %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid)))
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored <- d_cat %>%
group_by(subid, parent,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_scored_boot <- d_cat_scored %>%
group_by(target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_scored,
aes(x = target_num, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
scale_x_continuous(breaks = seq(0, 60, 12)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)")
ggplot(d_cat_scored,
aes(x = target_num, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root-transformed",
x = "age after square-root transformation (months)")
ggplot(d_cat_scored,
aes(x = target_ord, y = score, color = factor_names)) +
facet_grid(~ factor_names) +
geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_scored_boot, color = "black",
aes(y = mean, group = factor_names)) +
geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set2", guide = "none") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Age as an ordinal variable",
x = "age (ordinal)")
# contrasts(d_cat_scored$target_ord) <- contr.poly(13)
# contrasts(d_cat_scored$factor) <- contr.sum(4)
contrasts(d_cat$target_ord) <- contr.poly(13)
contrasts(d_cat$factor) <- contr.sum(4)
r1 <- lmer(response ~ factor * poly(target_num, 3) +
(1 + factor + target_num | subid),
d_cat %>%
mutate(target_num = scale(target_num, center = F)))
summary(r1)
Linear mixed model fit by REML ['lmerMod']
Formula: response ~ factor * poly(target_num, 3) + (1 + factor + target_num |
subid)
Data: d_cat %>% mutate(target_num = scale(target_num, center = F))
REML criterion at convergence: 536724.6
Scaled residuals:
Min 1Q Median 3Q Max
-5.5360 -0.4298 0.0225 0.5167 5.5499
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 200.49 14.159
factor1 156.13 12.495 0.65
factor2 145.83 12.076 -0.44 -0.59
factor3 69.75 8.352 -0.55 -0.34 -0.24
target_num 51.83 7.199 -0.64 -0.42 0.50 0.31
Residual 334.01 18.276
Number of obs: 61620, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) 68.3739 0.7634 89.569
factor1 6.9512 0.8216 8.461
factor2 -38.7255 0.7947 -48.728
factor3 25.1428 0.5573 45.116
poly(target_num, 3)1 3241.2664 88.4095 36.662
poly(target_num, 3)2 -1237.5812 18.2759 -67.717
poly(target_num, 3)3 532.0777 18.2759 29.114
factor1:poly(target_num, 3)1 -544.0136 31.6548 -17.186
factor2:poly(target_num, 3)1 2670.8981 31.6548 84.376
factor3:poly(target_num, 3)1 -2550.8843 31.6548 -80.584
factor1:poly(target_num, 3)2 32.5828 31.6548 1.029
factor2:poly(target_num, 3)2 38.2604 31.6548 1.209
factor3:poly(target_num, 3)2 881.6516 31.6548 27.852
factor1:poly(target_num, 3)3 48.0128 31.6548 1.517
factor2:poly(target_num, 3)3 -391.9779 31.6548 -12.383
factor3:poly(target_num, 3)3 -354.1489 31.6548 -11.188
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
This is weirdly hard to model… lots of alternative models didn’t converge. Also, those extremely high t-values make me a little suspicious, as do the identical standard errors for the last 9 terms. So, let’s take this with a bit of a grain of salt for now.
d_parent <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(d_demo %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid)))
Joining, by = "subid"
parent_counts <- d_demo %>%
distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
count(parent) %>%
mutate(proportion = n/sum(n))
knitr::kable(parent_counts, digits = 3)
| parent | n | proportion |
|---|---|---|
| No | 146 | 0.616 |
| Yes | 88 | 0.371 |
| NA | 3 | 0.013 |
d_demo %>%
count(Parent, ChildrenOldestAge_collapse)
Given these numbers, I think our best bet it just to look at parents vs. non-parents, and not try to separate out parents of young children (too few!).
scores_S2 <- efa_S2$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(d_parent %>% distinct(subid, parent)) %>%
left_join(d_cat %>% distinct(factor, factor_names)) %>%
mutate(factor = factor(factor))
Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector
# bootstrapped means and 95% CIs
d_parent_boot <- scores_S2 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5, # note: too close to dodge
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Exact age, untransformed",
color = "parent status: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent),
position = position_dodge(width = 0.3)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.3)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Exact age, square-root transformed",
color = "parent status: ",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_ord, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot, aes(y = mean, group = parent),
position = position_dodge(width = 0.5)) +
geom_pointrange(data = d_parent_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.5)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status",
subtitle = "Age as ordinal variable",
color = "parent status: ",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
# bootstrapped means and 95% CIs
d_cat_parent_scored_boot <- d_cat_scored %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_parent_scored_boot,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_parent_scored_boot,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformation",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_parent_scored_boot,
aes(x = target_ord, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot,
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Age as an ordinal variable",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
All of these visualizations suggest that, relative to non-parents (n = NA), parents (n = NA) tended to perceive greater abilities, in all domains except for bodily sensations. For both negative emotions and positive/social emotions, this difference is greatest in the mid-range of the target ages (12-24 months); for cognition & control, it’s greater at the older end of the target ages (e.g., 3-5 years). It doesn’t look like a huge effect, but it’s intruiging.
There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).
contrasts(scores_S2$parent) <- contr.treatment(2, base = 1) # baseline: NON-parents
contrasts(scores_S2$factor) <- contr.sum(4)
r2 <- lmer(score ~ target_num * factor * parent
+ (target_num + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r2)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * parent + (target_num + factor |
subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 22480.5
Scaled residuals:
Min 1Q Median 3Q Max
-9.3553 -0.4388 0.0524 0.5114 5.2974
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.3491 0.5908
target_num 0.0165 0.1284 -0.73
factor1 0.1947 0.4413 0.49 -0.40
factor2 0.2303 0.4799 -0.66 0.56 -0.40
factor3 0.2605 0.5104 0.24 0.00 -0.22 -0.61
Residual 0.2919 0.5403
Number of obs: 12168, groups: subid, 234
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.424399 0.049602 -8.556
target_num 0.279097 0.011334 24.624
factor1 0.078297 0.039260 1.994
factor2 -0.273518 0.042253 -6.473
factor3 0.253010 0.044634 5.669
parent2 0.070601 0.080884 0.873
target_num:factor1 -0.058596 0.006811 -8.603
target_num:factor2 0.190717 0.006811 28.001
target_num:factor3 -0.170985 0.006811 -25.104
target_num:parent2 0.003333 0.018482 0.180
factor1:parent2 0.005619 0.064020 0.088
factor2:parent2 -0.051089 0.068900 -0.741
factor3:parent2 0.031692 0.072783 0.435
target_num:factor1:parent2 -0.002939 0.011107 -0.265
target_num:factor2:parent2 0.053753 0.011107 4.840
target_num:factor3:parent2 -0.046204 0.011107 -4.160
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
The last three coefficients here are where the action is: According to this model, the difference between parents and non-parents in perceived developmental trajectories (target_num:parent2, which is not significant collapsing across factors) is exaggerated in the domain of cognition and control (target_num:factor2:parent2) and diminished in the domain of bodily sensations (target_num:factor3:parent2).
Let’s try adding polynomial effects:
r3 <- lmer(score ~ poly(target_num, 3) * factor * parent
+ (poly(target_num, 1) + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r3)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * parent + (poly(target_num,
1) + factor | subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 19808.5
Scaled residuals:
Min 1Q Median 3Q Max
-10.3042 -0.4652 0.0017 0.5305 5.8317
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.2278 0.4773
poly(target_num, 1) 519.8012 22.7991 -0.51
factor1 0.1982 0.4452 0.45 -0.39
factor2 0.2338 0.4835 -0.60 0.54 -0.40
factor3 0.2640 0.5138 0.30 0.00 -0.22 -0.61
Residual 0.2316 0.4812
Number of obs: 12168, groups: subid, 234
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.030561 0.039884 -0.766
poly(target_num, 3)1 48.855288 1.983984 24.625
poly(target_num, 3)2 -19.550418 0.613116 -31.887
poly(target_num, 3)3 8.448995 0.613116 13.780
factor1 -0.004388 0.038065 -0.115
factor2 -0.004395 0.041145 -0.107
factor3 0.011731 0.043587 0.269
parent2 0.075304 0.065038 1.158
poly(target_num, 3)1:factor1 -10.257028 1.061949 -9.659
poly(target_num, 3)2:factor1 3.413024 1.061949 3.214
poly(target_num, 3)3:factor1 -1.785482 1.061949 -1.681
poly(target_num, 3)1:factor2 33.384462 1.061949 31.437
poly(target_num, 3)2:factor2 3.384910 1.061949 3.187
poly(target_num, 3)3:factor2 -6.433547 1.061949 -6.058
poly(target_num, 3)1:factor3 -29.930494 1.061949 -28.184
poly(target_num, 3)2:factor3 8.594647 1.061949 8.093
poly(target_num, 3)3:factor3 -3.251241 1.061949 -3.062
poly(target_num, 3)1:parent2 0.583455 3.235228 0.180
poly(target_num, 3)2:parent2 -2.031708 0.999792 -2.032
poly(target_num, 3)3:parent2 1.883539 0.999792 1.884
factor1:parent2 0.001472 0.062072 0.024
factor2:parent2 0.024763 0.067094 0.369
factor3:parent2 -0.033507 0.071076 -0.471
poly(target_num, 3)1:factor1:parent2 -0.514473 1.731690 -0.297
poly(target_num, 3)2:factor1:parent2 -0.041808 1.731690 -0.024
poly(target_num, 3)3:factor1:parent2 1.099314 1.731690 0.635
poly(target_num, 3)1:factor2:parent2 9.409337 1.731690 5.434
poly(target_num, 3)2:factor2:parent2 -2.361622 1.731690 -1.364
poly(target_num, 3)3:factor2:parent2 -1.057673 1.731690 -0.611
poly(target_num, 3)1:factor3:parent2 -8.087893 1.731690 -4.671
poly(target_num, 3)2:factor3:parent2 6.669873 1.731690 3.852
poly(target_num, 3)3:factor3:parent2 -4.830063 1.731690 -2.789
Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
A lot to sort through here, but here are some observations:
poly(target_num, 3)2:parent2 and poly(target_num, 3)3:parent2)poly(target_num, 3)1:factor2:parent2)poly(target_num, 3)1:factor3:parent2). I’m not sure how to interpret the differences in non-linearity here.You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.
parent_counts_S1 <- demo_S1 %>%
distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
count(parent) %>%
mutate(proportion = n/sum(n))
knitr::kable(parent_counts_S1, digits = 3)
| parent | n | proportion |
|---|---|---|
| No | 147 | 0.653 |
| Yes | 77 | 0.342 |
| NA | 1 | 0.004 |
scores_S1 <- efa_S1$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(demo_S1 %>% distinct(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid))) %>%
mutate(factor_names = recode_factor(factor,
"MR4" = "Negative emotions (S1 F4)",
"MR1" = "Cognition & control (S1 F1)",
"MR2" = "Bodily sensations (S1 F2)",
"MR3" = "Positive/social emotions (S1 F3)"),
factor = factor(factor))
Joining, by = "subid"
# bootstrapped means and 95% CIs
d_parent_boot_S1 <- scores_S1 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S1 %>%
filter(!is.na(parent)) %>%
ggplot(aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_parent_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, group = parent)) +
geom_pointrange(data = d_parent_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = parent), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by parent status (STUDY 1)",
subtitle = "Exact age, untransformed",
color = "parent status: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
I won’t bother to re-plot with different treatments of age, since there are only three target ages here.
# make new dataframe
d_cat_S1 <- d_all_S1 %>%
rename(subid_target = X) %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(demo_S1 %>%
select(ResponseId, Parent) %>%
rename(subid = ResponseId, parent = Parent) %>%
mutate(subid = as.character(subid))) %>%
filter(!is.na(factor)) # drop extra items from S1 not in S2
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored_S1 <- d_cat_S1 %>%
group_by(subid, parent,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_parent_scored_boot_S1 <- d_cat_scored_S1 %>%
filter(!is.na(parent)) %>%
group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_parent_scored_boot_S1,
aes(x = target_num, y = score, color = parent)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_parent_scored_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, group = parent)) +
geom_pointrange(data = d_cat_parent_scored_boot_S1,
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
Looking quickly at these two visualizations of Study 1, I’d say that there’s a hint of the same patterns in the domains of negative emotions and positive social emotions (parents attributed more than non-parents, especially at 9 months), and in cognition/control (parents attributed more than non-parents at 5 years). So, generally consistent?
I won’t run regression analyses just now.
d_gender <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(d_demo %>%
select(ResponseId, GenderSex) %>%
rename(subid = ResponseId, gender = GenderSex) %>%
mutate(subid = as.character(subid)))
Joining, by = "subid"
gender_counts <- d_demo %>%
distinct(ResponseId, GenderSex) %>%
rename(subid = ResponseId, gender = GenderSex) %>%
count(gender) %>%
mutate(proportion = n/sum(n))
knitr::kable(gender_counts, digits = 3)
| gender | n | proportion |
|---|---|---|
| Female | 123 | 0.519 |
| Male | 114 | 0.481 |
scores_S2 <- efa_S2$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(d_gender %>% distinct(subid, gender)) %>%
left_join(d_cat %>% distinct(factor, factor_names)) %>%
mutate(factor = factor(factor), gender = factor(gender))
Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector
# bootstrapped means and 95% CIs
d_gender_boot <- scores_S2 %>%
filter(!is.na(gender)) %>%
group_by(gender, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S2 %>%
filter(!is.na(gender)) %>%
ggplot(aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_gender_boot, aes(y = mean, group = gender)) +
geom_pointrange(data = d_gender_boot, fatten = 1.5, # note: too close to dodge
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = gender), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant gender",
subtitle = "Exact age, untransformed",
color = "gender: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(gender)) %>%
ggplot(aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_gender_boot, aes(y = mean, group = gender),
position = position_dodge(width = 0.3)) +
geom_pointrange(data = d_gender_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.3)) +
# geom_smooth(aes(group = gender), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant gender",
subtitle = "Exact age, square-root transformed",
color = "participant gender: ",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(gender)) %>%
ggplot(aes(x = target_ord, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_gender_boot, aes(y = mean, group = gender),
position = position_dodge(width = 0.5)) +
geom_pointrange(data = d_gender_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.5)) +
# geom_smooth(aes(group = gender), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant gender",
subtitle = "Age as ordinal variable",
color = "participant gender: ",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
# bootstrapped means and 95% CIs
d_cat_gender_scored_boot <- d_cat_scored %>%
left_join(d_demo %>% distinct(ResponseId, GenderSex) %>%
mutate(ResponseId = as.character(ResponseId)) %>%
rename(subid = ResponseId, gender = GenderSex)) %>%
filter(!is.na(gender)) %>%
group_by(gender, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
Joining, by = "subid"
ggplot(d_cat_gender_scored_boot,
aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot,
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_gender_scored_boot,
aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot,
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformation",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_gender_scored_boot,
aes(x = target_ord, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot,
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Age as an ordinal variable",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
All of these visualizations suggest that, relative to women, men tended to perceive fewer abilities in the negative emotions domain (across the age range) and in the bodily sensations domain (especially for very young infants), and perhaps in the positive/social emotions domain (for toddlers). In contrast, if anythign men seemed to perceive more cognition & control abilities in the middle age ranges (~ 1 year). None of these are huge effects, but they are intriguing.
There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).
contrasts(scores_S2$factor)
[,1] [,2] [,3]
MR1 1 0 0
MR2 0 1 0
MR3 0 0 1
MR4 -1 -1 -1
r2 <- lmer(score ~ target_num * factor * gender
+ (target_num + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r2)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * gender + (target_num + factor |
subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 22684.2
Scaled residuals:
Min 1Q Median 3Q Max
-9.3855 -0.4361 0.0520 0.5098 5.3101
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.34624 0.5884
target_num 0.01648 0.1284 -0.73
factor1 0.18436 0.4294 0.49 -0.41
factor2 0.22582 0.4752 -0.65 0.57 -0.39
factor3 0.25876 0.5087 0.24 0.00 -0.23 -0.61
Residual 0.29011 0.5386
Number of obs: 12324, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.352499 0.053820 -6.550
target_num 0.279569 0.012339 22.658
factor1 0.198080 0.041761 4.743
factor2 -0.398918 0.045618 -8.745
factor3 0.304166 0.048464 6.276
gender2 -0.085230 0.077601 -1.098
target_num:factor1 -0.073166 0.007398 -9.890
target_num:factor2 0.232390 0.007398 31.413
target_num:factor3 -0.205753 0.007398 -27.813
target_num:gender2 -0.001484 0.017791 -0.083
factor1:gender2 -0.236684 0.060213 -3.931
factor2:gender2 0.211000 0.065775 3.208
factor3:gender2 -0.080956 0.069878 -1.159
target_num:factor1:gender2 0.028011 0.010667 2.626
target_num:factor2:gender2 -0.044941 0.010667 -4.213
target_num:factor3:gender2 0.037002 0.010667 3.469
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
The last six coefficients here are where the action is: According to this model, the difference between male and female in overall mental capacity attributions (gender2, which is not significant collapsing across factors and targets) is exaggerated in the domain of negative emotions (factor1:gender2), diminished in the domain of cognition and control (factor2:gender2), and neither diminshed nor exaggerated in the domain of bodily sensations (factor3:gender2). Meanwhile, the difference between male and female in perceived developmental trajectories (target_num:gender2, which is not significant collapsing across factors) is (likewise) exaggerated in the domain of negative emotions (target_num:factor1:gender2), diminished in the domain of cognition and control (target_num:factor2:gender2) and exaggerated in the domain of bodily sensations (target_num:factor3:gender2).
Let’s try adding polynomial effects:
r3 <- lmer(score ~ poly(target_num, 3) * factor * gender
+ (poly(target_num, 1) + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r3)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * gender + (poly(target_num,
1) + factor | subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 20032.2
Scaled residuals:
Min 1Q Median 3Q Max
-10.2922 -0.4698 0.0025 0.5288 5.8429
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.2251 0.4745
poly(target_num, 1) 519.0181 22.7820 -0.51
factor1 0.1878 0.4333 0.45 -0.40
factor2 0.2292 0.4788 -0.59 0.56 -0.39
factor3 0.2622 0.5120 0.29 0.00 -0.23 -0.61
Residual 0.2312 0.4809
Number of obs: 12324, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.04200 0.04320 0.972
poly(target_num, 3)1 48.93792 2.15991 22.657
poly(target_num, 3)2 -19.76380 0.66749 -29.609
poly(target_num, 3)3 8.89729 0.66749 13.329
factor1 0.09484 0.04043 2.345
factor2 -0.07099 0.04441 -1.599
factor3 0.01383 0.04733 0.292
gender2 -0.08732 0.06229 -1.402
poly(target_num, 3)1:factor1 -12.80752 1.15613 -11.078
poly(target_num, 3)2:factor1 2.82583 1.15613 2.444
poly(target_num, 3)3:factor1 -0.97761 1.15613 -0.846
poly(target_num, 3)1:factor2 40.67937 1.15613 35.186
poly(target_num, 3)2:factor2 3.21977 1.15613 2.785
poly(target_num, 3)3:factor2 -7.38714 1.15613 -6.390
poly(target_num, 3)1:factor3 -36.01663 1.15613 -31.153
poly(target_num, 3)2:factor3 13.88965 1.15613 12.014
poly(target_num, 3)3:factor3 -7.53568 1.15613 -6.518
poly(target_num, 3)1:gender2 -0.25972 3.11428 -0.083
poly(target_num, 3)2:gender2 -0.60651 0.96243 -0.630
poly(target_num, 3)3:gender2 0.27380 0.96243 0.284
factor1:gender2 -0.19716 0.05830 -3.382
factor2:gender2 0.14758 0.06403 2.305
factor3:gender2 -0.02874 0.06824 -0.421
poly(target_num, 3)1:factor1:gender2 4.90328 1.66698 2.941
poly(target_num, 3)2:factor1:gender2 0.95398 1.66698 0.572
poly(target_num, 3)3:factor1:gender2 -0.58404 1.66698 -0.350
poly(target_num, 3)1:factor2:gender2 -7.86688 1.66698 -4.719
poly(target_num, 3)2:factor2:gender2 -1.80751 1.66698 -1.084
poly(target_num, 3)3:factor2:gender2 1.33095 1.66698 0.798
poly(target_num, 3)1:factor3:gender2 6.47711 1.66698 3.886
poly(target_num, 3)2:factor3:gender2 -5.73172 1.66698 -3.438
poly(target_num, 3)3:factor3:gender2 5.11585 1.66698 3.069
Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
A lot to sort through here, but here are some observations:
poly(target_num, 2)1:factor2:gender2) and its cubic (poly(target_num, 3)1:factor2:gender2) components. But I’m having a hard time interpreting the coefficients at this point.You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.
gender_counts_S1 <- demo_S1 %>%
distinct(ResponseId, GenderSex) %>%
rename(subid = ResponseId, gender = GenderSex) %>%
count(gender) %>%
mutate(proportion = n/sum(n))
knitr::kable(gender_counts_S1, digits = 3)
| gender | n | proportion |
|---|---|---|
| Female | 96 | 0.427 |
| Male | 127 | 0.564 |
| Other (please describe) | 1 | 0.004 |
| Prefer not to say | 1 | 0.004 |
scores_S1 <- efa_S1$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(demo_S1 %>% distinct(ResponseId, GenderSex) %>%
rename(subid = ResponseId, gender = GenderSex) %>%
mutate(subid = as.character(subid))) %>%
mutate(factor_names = recode_factor(factor,
"MR4" = "Negative emotions (S1 F4)",
"MR1" = "Cognition & control (S1 F1)",
"MR2" = "Bodily sensations (S1 F2)",
"MR3" = "Positive/social emotions (S1 F3)"),
factor = factor(factor))
Joining, by = "subid"
# bootstrapped means and 95% CIs
d_gender_boot_S1 <- scores_S1 %>%
filter(!is.na(gender)) %>%
group_by(gender, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S1 %>%
filter(gender %in% c("Female", "Male")) %>%
ggplot(aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_gender_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 5),
aes(y = mean, group = gender)) +
geom_pointrange(data = d_gender_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = gender), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant gender (STUDY 1)",
subtitle = "Exact age, untransformed",
color = "participant gender: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
I won’t bother to re-plot with different treatments of age, since there are only three target ages here.
# make new dataframe
d_cat_S1 <- d_all_S1 %>%
rename(subid_target = X) %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(demo_S1 %>%
select(ResponseId, GenderSex) %>%
rename(subid = ResponseId, gender = GenderSex) %>%
mutate(subid = as.character(subid))) %>%
filter(!is.na(factor)) # drop extra items from S1 not in S2
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored_S1 <- d_cat_S1 %>%
group_by(subid, gender,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_gender_scored_boot_S1 <- d_cat_scored_S1 %>%
filter(!is.na(gender)) %>%
group_by(gender, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 5),
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
aes(x = target_num, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 1),
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 1),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformed",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
aes(x = target_ord, y = score, color = gender)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 0.5),
aes(y = mean, group = gender)) +
geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
filter(gender %in% c("Female", "Male")),
position = position_dodge(width = 0.5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
# scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "age (ordinal)",
x = "target",
y = "score (NOTE: scales differ)")
Looking quickly at these two visualizations of Study 1, I’d say that they are similar to the visualizations for Study 2, both in the negative emotions and bodily domains. There a hint of more of a difference here in the positive & social emotions domain. So, generally consistent?
I won’t run regression analyses just now.
d_edu <- d_all %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(d_demo %>%
select(ResponseId, Education) %>%
mutate(edu = case_when(as.numeric(Education) < 8 ~ "less than a BA",
as.numeric(Education) >= 8 ~ "BA or more"),
edu = factor(edu,
levels = c("less than a BA",
"BA or more"))) %>%
rename(subid = ResponseId) %>%
mutate(subid = as.character(subid)))
Joining, by = "subid"
edu_counts <- d_demo %>%
distinct(ResponseId, Education) %>%
rename(subid = ResponseId, edu = Education) %>%
count(edu) %>%
mutate(proportion = n/sum(n))
knitr::kable(edu_counts, digits = 3)
| edu | n | proportion |
|---|---|---|
| Some high school, no diploma | 2 | 0.008 |
| High school graduate, diploma or equivalent (including GED) | 31 | 0.131 |
| Some college credit, no degree | 61 | 0.257 |
| Trade school, technical school, or vocational school | 11 | 0.046 |
| Associate's degree (for example, AA, AS) | 26 | 0.110 |
| Bachelor's degree (for example, BA, BS) | 80 | 0.338 |
| Master's degree (for example, MA, MS) | 25 | 0.105 |
| Doctor or professional degree (for example, PhD, JD, MD, MBA) | 1 | 0.004 |
scores_S2 <- efa_S2$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(d_edu %>% distinct(subid, edu)) %>%
left_join(d_cat %>% distinct(factor, factor_names)) %>%
mutate(factor = factor(factor), edu = factor(edu))
Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector
# bootstrapped means and 95% CIs
d_edu_boot <- scores_S2 %>%
filter(!is.na(edu)) %>%
group_by(edu, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S2 %>%
filter(!is.na(edu)) %>%
ggplot(aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_edu_boot, aes(y = mean, group = edu)) +
geom_pointrange(data = d_edu_boot, fatten = 1.5, # note: too close to dodge
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = edu), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant edu",
subtitle = "Exact age, untransformed",
color = "edu: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(edu)) %>%
ggplot(aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_edu_boot, aes(y = mean, group = edu),
position = position_dodge(width = 0.3)) +
geom_pointrange(data = d_edu_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.3)) +
# geom_smooth(aes(group = edu), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant edu",
subtitle = "Exact age, square-root transformed",
color = "participant edu: ",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
scores_S2 %>%
filter(!is.na(edu)) %>%
ggplot(aes(x = target_ord, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_edu_boot, aes(y = mean, group = edu),
position = position_dodge(width = 0.5)) +
geom_pointrange(data = d_edu_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper),
position = position_dodge(width = 0.5)) +
# geom_smooth(aes(group = edu), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant edu",
subtitle = "Age as ordinal variable",
color = "participant edu: ",
x = "age (ordinal)",
y = "score (NOTE: scales differ)")
# bootstrapped means and 95% CIs
d_cat_edu_scored_boot <- d_cat_scored %>%
left_join(d_demo %>% distinct(ResponseId, Education) %>%
mutate(ResponseId = as.character(ResponseId),
edu = case_when(as.numeric(Education) < 8 ~ "less than a BA",
as.numeric(Education) >= 8 ~ "BA or more"),
edu = factor(edu,
levels = c("less than a BA",
"BA or more"))) %>%
rename(subid = ResponseId)) %>%
filter(!is.na(edu)) %>%
group_by(edu, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
Joining, by = "subid"
ggplot(d_cat_edu_scored_boot,
aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot,
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_edu_scored_boot,
aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot,
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformation",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_edu_scored_boot,
aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4) +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot,
position = position_dodge(width = 0.1),
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
position = position_dodge(width = 0.1),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformation",
x = "age after square-root transformation (months)",
y = "score (0-100)")
All of these visualizations suggest that, relative to less educated participants, more educated participants tended to perceive slightly fewer abilities in the negative emotions domain, but were otherwise generally very similar.
There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).
contrasts(scores_S2$edu) <- contr.treatment(2, base = 1) # baseline: less educated
contrasts(scores_S2$factor) <- contr.sum(4)
r2 <- lmer(score ~ target_num * factor * edu
+ (target_num + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r2)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * edu + (target_num + factor | subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 22726.8
Scaled residuals:
Min 1Q Median 3Q Max
-9.4244 -0.4335 0.0513 0.5090 5.3716
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.34804 0.5899
target_num 0.01646 0.1283 -0.73
factor1 0.19408 0.4406 0.49 -0.40
factor2 0.23103 0.4807 -0.66 0.56 -0.41
factor3 0.25862 0.5085 0.24 0.00 -0.21 -0.61
Residual 0.29081 0.5393
Number of obs: 12324, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.391254 0.052284 -7.483
target_num 0.282154 0.011952 23.607
factor1 0.081491 0.041379 1.969
factor2 -0.307750 0.044657 -6.891
factor3 0.290749 0.046956 6.192
edu2 -0.005014 0.078179 -0.064
target_num:factor1 -0.054535 0.007177 -7.598
target_num:factor2 0.208371 0.007177 29.033
target_num:factor3 -0.194921 0.007177 -27.159
target_num:edu2 -0.007374 0.017871 -0.413
factor1:edu2 0.006130 0.061873 0.099
factor2:edu2 0.023088 0.066774 0.346
factor3:edu2 -0.057068 0.070211 -0.813
target_num:factor1:edu2 -0.011532 0.010732 -1.075
target_num:factor2:edu2 0.005372 0.010732 0.501
target_num:factor3:edu2 0.015575 0.010732 1.451
Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
Nada.
Let’s try adding polynomial effects:
r3 <- lmer(score ~ poly(target_num, 3) * factor * edu
+ (poly(target_num, 1) + factor | subid),
scores_S2 %>%
mutate(target_num = target_num/12))
summary(r3)
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * edu + (poly(target_num,
1) + factor | subid)
Data: scores_S2 %>% mutate(target_num = target_num/12)
REML criterion at convergence: 20092.9
Scaled residuals:
Min 1Q Median 3Q Max
-10.3685 -0.4600 0.0056 0.5311 5.9055
Random effects:
Groups Name Variance Std.Dev. Corr
subid (Intercept) 0.2270 0.4764
poly(target_num, 1) 518.4274 22.7690 -0.51
factor1 0.1975 0.4444 0.45 -0.39
factor2 0.2344 0.4842 -0.60 0.55 -0.41
factor3 0.2620 0.5119 0.29 0.00 -0.22 -0.61
Residual 0.2322 0.4818
Number of obs: 12324, groups: subid, 237
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.006897 0.042031 0.164
poly(target_num, 3)1 49.390345 2.092247 23.606
poly(target_num, 3)2 -21.317213 0.648102 -32.892
poly(target_num, 3)3 10.204882 0.648102 15.746
factor1 0.004537 0.040120 0.113
factor2 -0.013717 0.043494 -0.315
factor3 0.015694 0.045850 0.342
edu2 -0.015421 0.062848 -0.245
poly(target_num, 3)1:factor1 -9.546136 1.122546 -8.504
poly(target_num, 3)2:factor1 2.830921 1.122546 2.522
poly(target_num, 3)3:factor1 -0.776043 1.122546 -0.691
poly(target_num, 3)1:factor2 36.474744 1.122546 32.493
poly(target_num, 3)2:factor2 3.825457 1.122546 3.408
poly(target_num, 3)3:factor2 -7.834359 1.122546 -6.979
poly(target_num, 3)1:factor3 -34.120438 1.122546 -30.396
poly(target_num, 3)2:factor3 11.638223 1.122546 10.368
poly(target_num, 3)3:factor3 -5.846249 1.122546 -5.208
poly(target_num, 3)1:edu2 -1.290868 3.128486 -0.413
poly(target_num, 3)2:edu2 2.820908 0.969091 2.911
poly(target_num, 3)3:edu2 -2.629115 0.969091 -2.713
factor1:edu2 -0.010143 0.059991 -0.169
factor2:edu2 0.030668 0.065035 0.472
factor3:edu2 -0.035090 0.068559 -0.512
poly(target_num, 3)1:factor1:edu2 -2.018615 1.678516 -1.203
poly(target_num, 3)2:factor1:edu2 1.014589 1.678516 0.604
poly(target_num, 3)3:factor1:edu2 -1.078782 1.678516 -0.643
poly(target_num, 3)1:factor2:edu2 0.940291 1.678516 0.560
poly(target_num, 3)2:factor2:edu2 -3.298146 1.678516 -1.965
poly(target_num, 3)3:factor2:edu2 2.431313 1.678516 1.448
poly(target_num, 3)1:factor3:edu2 2.726340 1.678516 1.624
poly(target_num, 3)2:factor3:edu2 -1.130454 1.678516 -0.673
poly(target_num, 3)3:factor3:edu2 1.724635 1.678516 1.027
Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
Again, nothing substantial - some indications of some education-related differences in perceptions of non-linear effects? But hard to make too much sense of that.
You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.
edu_counts_S1 <- demo_S1 %>%
distinct(ResponseId, Education) %>%
mutate(edu = case_when(
as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
rename(subid = ResponseId) %>%
count(edu) %>%
mutate(proportion = n/sum(n))
knitr::kable(edu_counts_S1, digits = 3)
| edu | n | proportion |
|---|---|---|
| less than a BA | 132 | 0.587 |
| BA or more | 93 | 0.413 |
scores_S1 <- efa_S1$scores %>%
data.frame() %>%
rownames_to_column("subid_target") %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(factor, score, -c(subid, starts_with("target"))) %>%
left_join(demo_S1 %>% distinct(ResponseId, Education) %>%
mutate(edu = case_when(
as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
rename(subid = ResponseId) %>%
mutate(subid = as.character(subid))) %>%
mutate(factor_names = recode_factor(factor,
"MR4" = "Negative emotions (S1 F4)",
"MR1" = "Cognition & control (S1 F1)",
"MR2" = "Bodily sensations (S1 F2)",
"MR3" = "Positive/social emotions (S1 F3)"),
factor = factor(factor))
Joining, by = "subid"
# bootstrapped means and 95% CIs
d_edu_boot_S1 <- scores_S1 %>%
filter(!is.na(edu)) %>%
group_by(edu, target, target_num, target_ord, factor_names) %>%
multi_boot_standard("score") %>%
ungroup() %>%
data.frame()
scores_S1 %>%
filter(!is.na(edu)) %>%
ggplot(aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.05) +
geom_line(data = d_edu_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 5),
aes(y = mean, group = edu)) +
geom_pointrange(data = d_edu_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = edu), #color = "black",
# method = "lm", formula = "y ~ poly(x, 3)") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
labs(title = "Developmental trajectories of factor scores, by participant education (STUDY 1)",
subtitle = "Exact age, untransformed",
color = "participant edu: ",
x = "target age (months)", y = "score (NOTE: scales differ)")
I won’t bother to re-plot with different treatments of age, since there are only three target ages here.
# make new dataframe
d_cat_S1 <- d_all_S1 %>%
rename(subid_target = X) %>%
mutate(subid = gsub("_.*$", "", subid_target),
target = gsub("^.*_", "", subid_target)) %>%
select(-subid_target) %>%
mutate(target_num = recode(target,
"target00mo" = 0,
"target0Xmo" = 4/30,
"target01mo" = 1,
"target02mo" = 2,
"target04mo" = 4,
"target06mo" = 6,
"target09mo" = 9,
"target12mo" = 12,
"target18mo" = 18,
"target24mo" = 24,
"target36mo" = 36,
"target48mo" = 48,
"target60mo" = 60),
target_ord = recode_factor(target,
"target00mo" = "newborns",
"target0Xmo" = "4-day-olds",
"target01mo" = "1-month-olds",
"target02mo" = "2-month-olds",
"target04mo" = "4-month-olds",
"target06mo" = "6-month-olds",
"target09mo" = "9-month-olds",
"target12mo" = "12-month-olds",
"target18mo" = "18-month-olds",
"target24mo" = "2-year-olds",
"target36mo" = "3-year-olds",
"target48mo" = "4-year-olds",
"target60mo" = "5-year-olds")) %>%
gather(capacity, response, -c(subid, starts_with("target"))) %>%
left_join(factors_S2) %>%
left_join(demo_S1 %>%
select(ResponseId, Education) %>%
mutate(edu = case_when(
as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
rename(subid = ResponseId) %>%
mutate(subid = as.character(subid))) %>%
filter(!is.na(factor)) # drop extra items from S1 not in S2
Joining, by = "capacity"
Joining, by = "subid"
And get an average “score” for each of these factors for each participant:
d_cat_scored_S1 <- d_cat_S1 %>%
group_by(subid, edu,
target, target_num, target_ord,
factor, factor_names) %>%
summarise(score = mean(response, na.rm = T)) %>%
ungroup()
# bootstrapped means and 95% CIs
d_cat_edu_scored_boot_S1 <- d_cat_scored_S1 %>%
filter(!is.na(edu)) %>%
group_by(edu, target, target_num, target_ord, factor, factor_names) %>%
multi_boot_standard("score") %>%
ungroup()
ggplot(d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 5),
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, untransformed",
x = "target age (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
aes(x = target_num, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 1),
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
filter(!is.na(edu)),
position = position_dodge(width = 1),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Paired") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "Exact age, square-root transformed",
x = "age after square-root transformation (months)",
y = "score (NOTE: scales differ)")
ggplot(d_cat_edu_scored_boot_S1 %>%
filter(edu %in% c("Female", "Male")),
aes(x = target_ord, y = score, color = edu)) +
facet_wrap(~ factor_names, ncol = 4, scales = "free") +
# geom_line(aes(group = subid), alpha = 0.15) +
geom_line(data = d_cat_edu_scored_boot_S1 %>%
filter(edu %in% c("Female", "Male")),
position = position_dodge(width = 0.5),
aes(y = mean, group = edu)) +
geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
filter(edu %in% c("Female", "Male")),
position = position_dodge(width = 0.5),
aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
# geom_smooth(aes(group = factor_names),
# method = "lm", formula = "y ~ poly(x, 3)",
# color = "black") +
scale_color_brewer(palette = "Set1") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "bottom") +
# scale_x_continuous(breaks = seq(0, 60, 12)) +
scale_y_continuous(breaks = seq(0, 100, 10)) +
labs(title = "Developmental trajectories of category 'summary scores'",
subtitle = "age (ordinal)",
x = "target",
y = "score (NOTE: scales differ)")
Looking quickly at these two visualizations of Study 1, I’d say that they are similar to the visualizations for Study 2 - perhaps a little more dramatic, especially in the negative emotions and bodily sensations domains.
I won’t run regression analyses just now.